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TIME  DEPENDENT  STUDIES  OF  REACTIVE  SHOCKS 
IN  THE  GAS  PHASE 

I . Introduction 

This  paper  preaents  reaulta  obtained  from  a detailed  numerical 
simulation  model  designed  to  aid  In  theoretical  and  experimental  Btudlas 
of  shock  and  detonation  phenomena.  The  calculations  describe  the 
non-equilibrium  chemical  kinetics,  reaction  products,  and  intermediates 
produced  in  reactive  gas  mixtures  ignited  by  the  propagation  of  a 
shock  front.  The  model  is  based  on  a rather  general,  compressible, 
one- dimensional  reactive  flow  model  designed  to  treat  multi-fluid 
problems.  Since  the  fluid  dynamics  and  the  chemical  kinetics  are  time- 
dependent  and  self-consistent,  we  can  observe  the  effects  of  the  chemical 
energy  release  on  the  fluid  flow  during  highly  exothermic  reactions  and 
the  effects  of  temperature  and  density  inhomogeneities  in  the  reaction 
medium  on  the  resulting  chemical  composition. 

In  order  to  illustrate  one  type  of  physical  problem  which  the  reactive 
shock  model  can  address,  we  simulate  a shock  tube  in  which  a diaphragm  is 
burst  at  the  onset  of  the  calculation.  We  assume  that  initial  tempera- 
tures have  equilibrated  across  the  diaphragm  which  separates  regions  of 
differing  pressures.  The  calculational  method  discussed  below  has  first 
been  tested  against  the  results  obtained  from  a derivation  of  the  shock 
parameters  in  an  ideal  gas  with  no  conductivity,  diffusion,  or  chemical 
kinetics.  Wa  then  proceed  to  examine  a shock  propagating  in  an  H2-og 
medium.  The  detailed  kinetics  are  confined  to  this  relatively  simple 
combustion  process  since  it  has  bean  studied  extensively  using  shock  tube 
aa  wall  as  flame-structure  experiments.  Wa  look  in  detail  at  two  reactive 
flow  cases:  one  In  which  the  chemical  energy  is  released  quickly  enough 

to  altar  the  properties  of  the  shook,  and  tht  second  in  which  the  shock 
Note:  Manuicript  submitted  October  4,  1978. 
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Is  reflected  from  a rigid  wall.  We  look  behind  the  shock  as  It  reflects 
off  the  end  wall  and  observe  the  increased  temperature  and  reaction  products. 

In  the  sections  that  follow,  the  numerical  algorithms  which  are  the 
foundation  of  the  fluid  transport  and  chemical  Integration  are  described, 
followed  by  a description  of  the  chemical  reaction  scheme  and  the  physical 
processes  Included  in  the  model.  In  Section  IV  several  calculations  are 
described.  The  initial  evolution  of  a shock  propagating  in  an  mixture 

with  and  without  chemical  reactions  is  compared  to  calculations  performed  in 
gas  mixtures  where  fluid  dynamic  motion  is  suppressed.  The  later  evolution  of 
the  reactive  shock  is  alao  presented  as  the  shock  reflects  from  a rigid  wall. 

Before  proceeding,  however,  it  is  useful  to  put  this  modal  and  the  sample 
calculation  in  perspective  with  other  types  of  calculations  used  for  shock  and 
detonation  studies.  First,  this  la  a time-dependent  model,  but  not  a time- 
dependent  kinetics  or  hydrodynamic  calculation  alone;  special  emphasis  hss 
been  placed  on  the  preclae  coupling  of  these  interactions.  Thus  the  model  is 
used  to  study  transient  phenomena,  the  approach  to  steady  state,  the  way  in 
which  a perturbation  on  the  system  alters  the  steady  state,  the  effects  of 
non-ideality  of  gases,  and  the  way  in  which  the  release  of  chemical  energy 
altera  flow  patterns.  Because  the  algorithm  used  to  solve  the  hydrodynamics 
has  minimal  numerical  diffusion  and  no  artificial  viscosity,  a gradient  as 
steep  as  a shock  front  can  be  resolved  in  one  or  two  cells.  Thus  tha  model 
could,  for  example,  be  used  to  perform  detailed,  time-dependent  simulations 
of  how  the  shape  of  the  ahock  front  can  ba  altered  by  endothermic, 


exothermic,  or  relaxation  processes.  In  summary,  the  major  advance  we  are 
describing  is  the  application  of  extremely  accurate,  efficient  state-of- 
the-art  numerical  techniques  to  the  study  of  the  interplay  between  funda- 
mental flow  and  chemical  processes. 
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IX.  Numerical  Modal 


We  solve  the  tims-depeodent  equetlona  for  conservation  of  mess, 
momentum,  and  energy! 
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where  p,  pv,  E and  P are  the  total  mass,  momentum,  and  energy 
densities  and  pressure,  respectively.  The  {p^}  and  { > are  the  maBa 
density  and  the  diffusion  velocities  of  the  Individual  species.  The 
quantities  T1  and  x represent  the  viscosity  and  the  thermal  conductivity 
of  the  gee  mixture  et  specified  {Pl},  and  temperature,  T.  The  and 
{Lt]  refer  to  chemical  production  end  lose  processes  for  species  1,  and 
the  lest  term  in  Eq.  (4)  represents  the  change  In  energy  due  to  chemlcel 
reactions  which  Is  added  to  the.  fluid  dynamic  energy  density.  This  term 
only  enters  the  total  energy  equation  Implicitly  and  is  tharafora  enclosed 
In  brackets.  The  model  we  have  used  Is  one -dimensional,  end  the  spatial 
gradients  era  evaluated  only  in  the  r-dlrectlon,  where  r is  a generalised 
position  coordinate.  Currently  we  have  assumed  Cartesian  coordinates  for 
our  shock  tube  simulations  but  the  model  allow*  for  other  coordinate  systems. 
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There  are  a number  of  inherent  numerical  difficulties  we  face  in 
detailed  modeling  of  the  transport  of  a shock  or  a flame  front  in  a 
reacting  medium.  One  problem  arises  as  the  result  of  the  widely 
different  time  scales  characteristic  of  fluid  and  chemical  processes, 
Another  obstacle  to  numerical  simulation  Is  that  conventional  numerical 
methods  are  uneble  to  resolve  accurately  the  steep  gradients  in  pres- 
sure, density,  and  temperature  characteristic  of  these  phenomena. 
Finally,  we  face  the  problem  of  a time-consuming  chemical  rate  inte- 
gration which,  because  of  the  many  species  and  reactions  required  to 
model  a realistic  system,  can  Increase  calculation  time  by  orders  of 
magnitude  over  a pure  fluid  calculation. 

We  approach  tha  time-scale  problem  by  employing  the  fractlonal- 
atep  method  commonly  known  as  time -step  splitting.  For  example,  the 
apecies  density  equations  are  separated  Into  spatial  derivative  terms 
and  chemical  production  and  losa  terms.  The  resultant  sat  of  partial 
differential  equations  are  solved  separately  from  the  set  of  nonlinear 
ordinary  differential  equations  at  each  tima  step.  The  solution  is  a 
result  of  coupling  the  fractional  steps,  accomplished  by  providing  the 
initial  conditions  for  ons  step  with  the  results  of  the  solution  from 
the  other  step.  This  technique  is  advantageous  since  very 
efficient  and  modular  computer  codes  can  be  produced  and  the 
various  parts  of  an  aquation  may  ba  solved  with  an  appropriate  time- 
step.  Accuracy  Is  testad  by  dacreasing  the  time- step  and  by  increasing 
the  resolution,  and  then  checking  for  variations  in  the  answers. 
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The  fluid  dynamic  equations  expreaging  the  conservation  of  mass, 

momentum,  and  energy  are  solved  by  the  Flux-Corrected  Transport  (FCT) 

12  3 

method,  ’ * as  are  the  continuity  equations  for  the  densities  of  the 
individual  molecular  species.  FCT  is  designed  to  reproduce  steep 
gradients  particularly  well,  a situation  which  occurs  fluid  dynamically 
in  shocks  and  chemically  in  combustion  fronts  and  flames.  Furthermore, 
because  negative  densities  cannot  occur  using  the  method,  many  types 
of  numerical  instabilities  which  plague  reactive  flow  calculations 
can  be  avoided  without  excess  non-physical  numerical  diffusion, 

The  particular  FCT  algorithm  chosen  is  a simple  variant  of  the 
ETBFCT  program , Instead  of  following  the  cell  center  locations  on 
the  generalized  one-dimensional  grid  and  averaging  these  to  define  cell 
boundaries,  we  follow  the  cell  boundaries  as  primal  quantities  and 
average  these  locations  to  define  cell  centers.  By  ascribing  physical 
reality  to  the  cell  boundaries,  numerical  diffusion  between  ran  be 
eliminated  and  l.agrangian  treatment  of  system  boundaries  it 
simplified. 

An  important  feature  of  these  FCT  algorithms  is  the  lack  of  an 
artificial  viscosity  needed  to  stabilize  ehocke.  The  flux  correction 
procedure  itself  ensures  that  the  shocks  are  one  or  two  zones  wide  and 
have  maximal  resolution.  Physically  correct  mass  diffusion,  viacoaity, 
and  thermal  conduction  terms  are  included  in  the  model  bye  separate 
diffusion  subprogram  which  takes  advantage  of  time-step  splitting.  The 
fluid  dynamics  time  integration  is  performed  by  an  explicit  two  step 
predictor-corrector  technique^.  The  usual  Courant-Friedrich-Lewy 


•tability  condition*  apply  to  the  hydrodynamics  time-step.  However  the 
chemistry  can  proceed  on  a faster  scale  provided  that  the  hydrodynamics 
does  not  respond  to  these  chemical  changes  on  this  faster  scale. 

Another  important  feature  of  the  FCT  algorithms  chosen  is  their 
ability  to  divorce  the  grid  motion  from  the  fluid  flow.  We  have  used 
this  freedom  to  include  an  efficient  adaptive  grlddlng  procedure  to 
automatically  follow  regions  such  as  the  shock  front,  where  enhanced 
resolution  is  required,  in  this  paper  the  adaptation  is  based  on  the 
shock.  The  region  innediately  around  the  shock  front  and  for  pre- 
determined distances  on  each  side  is  gridded  with  finely  spaced  cells. 
The  fine  spacing  transitions  smoothly  into  the  more  coarsely 
resolved  region.  As  the  shock  moves  along  the  length  of  the  tube,  the 
finely-spaced  region  moves  with  it  and  reflects  off  the  boundary  wall. 

We  have  found  that  the  condition  on  the  acceleration, 

I LS  I = maximum, 

1 Ar  1 £ 

is  adequate  for  locating  the  shock  front,  Here  P is  the  pressure,  r 
is  the  generalized  position  coordinate,  and  "p  is  the  average  maas 
density.  In  principle  there  can  be  any  number  of  finely- spaced  regions 
and  the  method  will  also  be  useful  for  obtaining  the  required  resolution 
at  a flame  front. 

The  rate  equations  for  the  species  number  densities  can  be  written 
as  a set  of  first-order,  coupled,  nonlinear, ordinary  differential 
equations  of  the  following  form  for  the  1 species: 
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where  Is  the  density,  Is  the  production  rate,  and  the 
characteristic  loss  time, 

b- 

The  "selected  asymptotic  method"  * used  to  solve  these  equations, 
determines  which  equations  satisfy  an  appropriately  chosen  stiffness 
criterion.  The  non-stiff  equations  are  solved  by  a second  order  Adams 
method  while  a very  stable  asymptotic  method  Is  applied  to  the  stiff 
equations.  This  method  la  particularly  suitable  for  the  situation 
which  frequently  occurs  in  the  near  equilibrium  state  when  the  production 
and  loss  rates  are  very  large  but  nearly  cancel.  Small  numerically 
produced  oscillations,  which  result  from  short  time  constants  in  the 
kinetics  and  would  produce  prohibitively  short  time  steps,  are  effec- 
tively reduced.  The  time  step  is  carefully  monitored  to  insure  the 
accuracy  and  convergence  of  the  method.  In  coupled  chemical  and  hydro- 
dynamic  models  the  accuracy  of  the  chemical  Integration  can  be  relaxed 
somewhat  due  to  the  presence  of  unavoidable  truncation  errors  in  the 
fluid  dynamics  and  transport  schemes  and  the  uncertainties  In  the  rate 
constants. 

The  selected  asymptotic  method  iu  very  efficient  for  reactive 
flow  problems  since  it  has  very  low  overhead  and  can  be  restarted  Inex- 
pensively. It  hea  been  programmed  to  make  efficient  uae  of  the  parallel 
processing  capability  of  vector  computers. 
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In  the  time-splitting  approximation,  wa  assume  that  the  fluid 
motion  Is  frozen  while  the  chemical  rate  equations  are  advanced  ualng 
their  own  appropriate  time-step.  At  the  end  of  a fluid  Integration 
step,  the  internal  energy,  £ , is  obtained  by  evaluating 


E- 
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(6) 


p Pv  . 

where  E,  pv,  and  p have  bean  determined.  Then  the  temperature,  which 
la  required  for  evaluation  of  chemical  reaction  rates  and  diffusion 
coefficients,  is  found  by  iterating  an  equation  of  the  form 


€ * h (T,  {ni} ) - N k T 


<7> 


where  C la  in  ergs/cm  , h ia  the  enthalpy,  [n^}  and  N are  the  special 
and  total,  number  density,  respectively,  and  k 1b  Boltzmann's  constant. 
Thua  during  a chemical  integration  stap,  T la  found  after  ths  equation* 
are  aolved  for  [n^].  ’Die  internal  energy  of  the  system  1b  held 
constant  for  each  chemical  time-step.  We  emphasize  that  at  each  time- 
step  the  enthalpy  and  the  specific  heats  are  re-evaluated  using  the  new 
densities  and  temperature. 


The  diffusion  velocities  {V^}.  are  found  by  solving 
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subject  to  the  constraint 


I P.V  - 0.  (9) 

J J J 

The  {D..}  are  binary  diffusion  coefficients  and  will  be  discussed  In 

J K 

the  next  section.  The  {DTj}  are  the  thermal  diffusion  coefficients 
and  have  be*n  set  equal  to  zero  in  the  calculations  presented  below. 
Equations  (8)  and  (9)  form  a set  of  linear  coupled  equations  in  { } 
which  are  solved  by  an  iterative  technique  which  vill  be  described 
in  detail  in  a future  paper. 


III.  Chemical  Modal 

Tha  results  presented  below  are  derived  from  elmuletlona  of 
shocks  In  a mixture  of  hydrogen  and  oxygen  which  may  Include  an 
Inert  species  such  as  argon  or  helium.  This  relatively  simple 
combustion  system  Is  Ideal  for  our  benchmark  work  because  it  has 
been  studied  extensively  In  laboratory  shook  tubes  and  because 
a great  deal  of  effort  by  the  combustion  community  has  gone  Into 
determining  the  relevant  reaction  mechanisms  and  rates  over  a wide 
temperature  rango^,^,^,t?,1°.  For  combustion,  we  monitor  at 

each  point  in  the  system  the  evolution  in  time  of  the  temperature 
and  eight  reacting  species  densities,  H^,  0^,  H^o,  II,  0,  OH,  H0g,  and 
HgO^, and  the  density  of  the  inert  material. 

There  are  over  forty  reaction  paths  Included  here  which  are 
of  varying  lmportanca  at  aach  tamparatura  and  which  control  tha 
intermediate  and  final  denaltias.  These  era  summarized  In  Table  I. 

The  main  features  of  the  Induction  and  termination  parlode  are  fairly 
wall  understood.  Useful  for  our  purposes  Is  the  fact  that  in  tha 
feat  reaction  region  In  which  wa  are  interacted,  tha  major  termination 
patha  are  through  thermoleculer  raactlone  and  not  wall  reactions. 

Wa  ere  relatively  secure  in  our  treatment  of  termolecular  reaction 
retea,  which  ere  properties  of  tha  composition  of  the  systems;  the 
well  reactions  would  vary  with  the  well  costings  end  require  empirical 
modeling.  In  the  work  presented  here,  we  assume  there  is  one  rate 
coefficient  for  e termolecular  reaction,  regardless  of  whet  the  third 
body  is.  There  is  some  information  evtilable  on  the  efficiencies  of 
verloue  species,  and  this  can  be  Incorporated  in  future  calculations.1* 


n 


The  enthalpies  for  the  species  we  monitor  have  been  taken  from  the 
12 

JAMAF  tables  and  from  the  sixth  order  enthalpy  coefficients  of  Gordon 
and  McBride.1^  For  each  species  va  have  transformed  their  values  Into 
a seventh  order  polynomial  fit  which  goes  from  0°K  to  the  maximum 
temperature  given  In  the  JANAF  tables.  The  enthalpy  of  a species  "1" 
can  then  be  written 


i 5 

where  h is  In  units  of  erga/cm  and  T is  In  dsgrees  Kelvin.  The 
coefficients  a*  are  given  in  Table  II.  Our  motlvlatlon  for  using  one 
fit  over  a wide  temperature  range  Instead  of  table  look-ups  or  different 
sets  of  coefficients  for  the  temperature  range  below  and  above  1000°K 
lias  in  the  speed  and  ease  with  which  the  evaluation  of  expressions 
such  as  Eq.  (10)  can  be  vectorised  for  efficient  computation. 

The  transport  processes^  such  as  diffusion,  thermal  conduction, 
and  viscous  dissipation  are  generally  slow  compared  to  the  fluid 
flow  velocities  with  which  we  are  concerned  in  a shock  problem  and 
their  effects  on  the  reaction  cone  are  predictably  small.  We  have,  how- 
ever, Included  them  in  the  model  for  completeness  and  because  this  generic 
model  will  also  be  applied  to  tha  slower  propagation  of  a flam# 
front.  We  will  comment  briefly  In  the  sections  below  on  the  also  of 
the  diffusion  end  conductivity  offsets.  The  transport  coefficients  we 
have  choaen  ere  e flrat  approximation.  Whan  this  modal  Is  applied  to  e 
propagating  flame  front,  w*  will  try  to  incorporate  more  accurate  values 
of  the  quantities. 
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Accurate  binary  diffusion  coefficients  are  given  by  Marrero  and 
Mason, ^ These  have  been  used  where  available.  When  a coefficient 
was  not  tabulated,  we  have  exploited  the  similarity  between  diffusion 
of  species  close  in  molecular  weight  to  approximate  the  expression 
we  need.  Table  111  summarizes  the  binary  diffusion  coefficients  we 
have  used. 

The  thermal  conductivity  of  species  "i"  is  evaluated  in  terms  of 
the  viscosity,  n^,  the  ratio  of  specific  heats,  and  the  molecular 
weight,  m^: 


X 


i 


1 


<9  Y£“5) 


8.12  x 107 
(v1“l)m1 


.Hi 


(id 


The  units  of  XA  are  erg-cm"1-sec  This  expression  reduces  to 

the  usual  monatomic  gas  expression  for  the  appropriate  value  of  y^. 
Values  of  T| ^ , summarized  in  Table  IV,  are  admittedly  only  first 
approximations. 

The  conductivity  of  the  gas  mixture  X^,  has  been  evaluated  from 
an  approximate  expression  of  the  form: 


_1_ 


K y n® 


(12) 


where  the  expression  for  T|M  is 

(1.3) 

Mi 

raiN  Mi  Dij 

We  note  that  Eq.  (10)  for  XM  is  purely  empirical  and  based  on  a general - 

Im- 
itation of  the  tvo-specles  empirical  formula  , After  examining 

measured  values  of  T|M  for  many  mixtures  we  decided  this  sxpresslon  gives 

eccurets  enough  veluee  for  the  density  and  temperature  ranges  with 

which  we  ere  concerned  here. 


l 
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IV.  Discussion  of  Results 


We  have  performed  a number  of  simulations  of  shocks  in  an  Hg-Og 
medium  varying  the  initial  temperatures,  pressures  and  amounts  of 
diluent  on  the  low  density,  reacting  side  of  the  diaphragm.  In  this 
section  we  present  calculations  of  a shock  propagating  through  a 
reacting  gas,  initially  at  room  temperature,  with  no  diluent  present 
to  absorb  heat.  Within  the  constraints  of  a one-dimensional  model, 
these  calculations  show  how  the  energy  released  during  the  chemical 
reactions  affects  the  fluid  motion. 

Figure  1 summarizes  the  initial  conditions  before  the  diaphragm 
is  burst.  The  driver  gas  is  helium  at  10  atm  and  the  low  density 
region  is  a stoichiometric  mixture  of  and  0^  at  0.01  atm.  Small 
amounts  of  the  intermediate  chain  centers,  the  peroxides,  and  water 
have  been  included  on  the  reacting  side.  The  initial  densities  of 
these  minor  constituents  have  been  determined  in  a separate  calcula- 
tion by  allowing  a stoichiometric  mixture  of  and  0^  to  "equilibrate" 
at  fixed  temperature  (300°K)  for  about  an  hour  of  physical  time.  The 
chemical  reactions  described  by  the  rates  in  Table  I control  this 
process. 

Table  V gives  the  solution  to  an  ideal  diaphragm-shock  problem 
defined  in  Fig.  1 and  is  the  analytic  solution  to  the  Rankine- 
Hugoniot  equations.  These  results  are  strictly  valid  only  when  there 
are  no  chemical  reactions,  diffusion,  or  conductivity  in  the  fluid  and 
when  the  gas  constant  is  independent  of  temperature.  Therefore  Table  V 
can  only  be  used  as  a guide  to  understanding  the  detailed  simulation  pre- 
sented below  in  which  chemical  reactions,  molecular  diffusion,  thermal 
conductivity,  and  the  temperature  dependence  of  the  specific  heats  are 
explicitly  included. 
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At  the  onset  of  the  calculation,  the  diaphragm  la  removed  and  the 

driver  fluid  begins  to  move  Into  the  low  density  region.  The  five 

distinct  regions  described  by  Table  V are  not  clearly  differentiated  In 

-5 

our  calculations  with  0,1  cm  cones  until  ~ 1x10  sec.  Then  we  begin 
to  aee  the  clear  formation  of  a shock  front,  contact  discontinuity,  and 
expansion  fan.  By  about  8x10  sec  the  fluid  flow  has  settled  down  and 
looks  very  much  like  the  ideal  pattern.  Until  enough  time  has  elapsed 
so  that  there  are  substantial  changes  In  the  species  densities  due  to 
the  chemical  reactions,  the  fluid  motion  is  relatively  unaffected  by 
chemistry.  Detailed  comparisons  over  a wide  range  of  shock  parameters 
during  this  early  time  have  verified  the  accuracy  of  the  non-reacting 
fluid  model  as  compared  to  the  ideal  conditions  in  Table  V.1 

Figure  2a  shows  the  temperature  profile  after  ~ 1. 3^x10 ’ sac. 

At  thla  time  the  shock  front  is  located  at  kO  cm  down  the  tube.  The 
different  regions  are  clearly  separated  in  the  calculation  and  the 
energy  released  by  chemical  reactions  has  noticeable  effects  on  both 
the  shock  speed  and  the  shape  of  the  shocked  region.  The  dashed  line 
marks  the  corresponding  temperature  profile  for  the  ideal  non-reacting 
fluid  calculation.  An  increase  in  the  temperature  at  the  back  of  the 
shockad  region,  where  most  of  the  chemical  reactions  have  occurred,  Is 
clearly  visible.  The  calculation  also  shows  a broadening  of  the  region 
surrounding  the  contact  discontinuity  due  to  diffusion  of  the  helium 
driver  gae  into  the  reactive  mixture.  The  effect  of  thie  temperature 
decrease  is  apparent  in  Fige.  2b,  2c  and  3, 


Figure  2b  shows  the  total  densities  corresponding  to  the  temperatures 
shown  in  Fig.  2u.  As  one  might  expect,  the  high  temperature  regime 
corresponds  to  a low  density  regime  in  the  shocked  region.  The  total 
pressure  corresponding  to  Figs,  2a  and  2b  is  shown  ir.  Fig.  2c.  We  note 
that  both  temperature  and  total  density  can  vary  by  as  much  aa  a factor 
of  two  in  the  shocked  region,  whereas  the  pressure  variation  is  at  moat 
about  30%.  The  calculated  fluid  velocity  has  not  been  shown  since  it 
varies  at  most  by  10%  from  the  solutions  presented  in  Table  V. 

Figure  3 hIiows  the  number  density  of  and  0^,  the  chain  centers, 

and  the.  product  H^O  as  a function  of  location  along  the  length  of  the 
— 

tube  at  1.47  X 10  sec.  if  we  were  looking  at  a case  where  the  tempera- 
turo  were  uniform  across  the  high  temperature  region  (i.e.,  large  amounts 
of  diluent  present),  we  would  expect  the  H^O  peak  to  be  in  front  of  the 
OH  peak.  In  such  a ease  the  spatial  difference  between  the  peaks  would 
bo  analogous  to  the  temporal  separation  butween  peak  OH  formation  which 
precedes  the  maximum  product  formation  In  a homogeneous  reacting  gas 
mixture.  Looking  for  this  affect  in  the  highly  exothermic  case  we  are 

studying  can  be  misleading:  here  we  have  to  consider  both  kinetic  and 

fluid  dynamic,  effects.  The  high  OH  densities  indicate  that  combustion 
is  occurring  between  about  36  and  38  cm  and  this  is  in  agreement  with 
the  high  temperatures  shown  in  Fig.  2a.  The  unreacted  and  Oj,  near 
35  cm  is  in  the  cool  region  near  the  contact  discontinuity  into  which 
helium  has  diffused. 

We  have  used  the  same  numerical  model  in  a slightly  different 

mode  to  calculate  the  reactions  In  a stationary  premlxud  system  aver  a 

range  of  temperatures.  This  is  done  by  by-passing  the  transport  part 
of  the  calculation,  setting  the  density  across  the  grid  equal  to  that 
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behind  the  shock  front,  and  at  the  onset  oi  the  calculation  putting  a 
linear  temperature  ramp  across  the  grid.  With  no  fluid  motion  and  no 
diffusion  processes,  this  is  equivalent  to  performing  a large  number 
of  independent  tine-dependent  chemical  rate  calculations  simultaneously. 

The  difference  between  highly  exothermic  chemical  reactions 
occurring  behind  a moving  shock  front  and  reactions  in  a stationary 
premixed  homogeneous  gas  is  shown  in  Fig.  4.  Using  the  temperature 
ramp  method  described  above  for  the  range  1100-l‘ji00oK,  we  have  moni- 
tored the  gas  mixture  as  it  is  auddenly  brought  to  high  temperatures, 
raacta,  and  aelf-conaiatently  heats.  Figure  4a  shows  the  self-consis- 
tent temperature  Increase  with  time.  The  curves  represent 
stoichiometric  llg-0g  mixtures  with  density  equal  to  that  immediately 
behind  the  shock  front.  Initially  these  mixtures  are  suddenly  heated 
to  1100,  1200,  and  l'jOO°K  respectively. 

The  three  curves  in  Fig.  4b  show  the  increases  in  the  HgC  density 
corresponding  to  the  temperature  profiles  in  Fig.  4a.  After  an 
ignition  parlod,  the  bulk  of  the  water  is  formed  and  the  cuxvea  level 
off.  Alao  on  this  figure  we  have  plotted  the  maximum  water  density 
calculated  behind  the  moving  shock  front  at  a function  of  time. 

Note  that  the  temperature  behind  the  ideal  shock  is  about  1200°K,  and 
the  shock  results  fall  between  the  stationary  1200  and  1?00°K  curves 
with  respect  to  ignition  times,  but  below  these  curve*  in  maximum  dens- 
ity values.  This  latter  difference  is  the  natural  expansion  which 
occurs  whan  chamlcal  anargy  is  relsssad. 

Another  striking  difference  between  the  homogeneous,  prsmixad  end 
the  shock  calculations  is  the  difference  in  maximum  temperature 
attained.  Behind  the  shock  this  la  about  2300“2400°K,  wheraaa  the 
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calculations  in  Fig.  4 show  a maximum  of  about  3200°K.  This  difference 
in  thermal  energy  Is  that  energy  used  to  accelerate  the  shock  front. 

At  this  point  we  note  that  the  calculations  show  the  presence  of  a 
temperature  oscillation  in  the  shocked  region  even  when  there  are  no 
chemical  reactions  permitted.  At  1.34x10"*'  sec,  these  waves  have  an 
amplitude  of  about  25°K,  and  they  are  damped  somewhat  when  thermal 
conductivity  Is  allowed.  These  are  small,  finite  amplitude  sound 
waves,  a result  of  initial  transients.  As  the  shock  region  spreads  out 
and  time  increases,  the  amplitude  of  these  waves  also  decrease.  The 
oscillations  are  primarily  in  the  fluid  velocity  and  this  in  turn  Is 
amplified  in  the  internal  energy  and  the  pressure  calculation  through 
Eq.  (6). 

The  second  problem  we  shall  discuss  is  that  of  a shock  reflected 
off  of  a rigid  wall.  Many  shock  tube  experiments  are  designed  to  take 
advantage  of  the  relatively  uniform  regime  near  the  wall  behind  the 
reflected  shock,  The  temperature  and  density  are  further  elevated 
here  and  the  fluid  velocity  is  essentially  aero.  Thus  there  is  a 
period  of  relative  calm  that  can  be  used  to  watch  reactions  develop. 

In  order  to  model  the  region  behind  a reflected  shock,  we  have 
allowed  our  finely  grldded  region  to  move  to  the  wall  with  the  initial 
shock  front.  Then  the  grid  motion  is  restricted  so  that  it 
does  not  move  backwards.  Thus  we  retain  fine  resolution  near  the  wall. 

We  look  at  the  same  shock  as  described  above  and  set  the  wall  at  30  cm. 

For  this  problem  we  present  our  results  a little  differently.  We 
assume  that  there  is  s "hole"  1 cm  from  the  reflection  wall,  at  49  cm, 
and  monitor  what  passes  by,  Figure  ‘j  shows  how  the  temperature  at  this 
location  changes  as  a function  of  time.  When  the  shock  passes  49  cm,  the 
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temperature  Jumps  Co  about  2300°K  which  in  Chla  case  is  characteristic 

of  the  temperature  behind  the  shock  after  burning.  The  temperature 

stays  at  this  value  until  the  reflected  shock  passes  over  this  position 

again.  Then  the  temperature  increases  to  about  3500°K  and  the  pressures 

increase  by  a factor  of  about  three.  These  increases  are  consistent 

with  what  is  expected  from  a solution  of  the  Jung>  conditions  for  an 

17 

ideal  shock  reflected  from  a wall 

The  calculated  fluid  velocity  behind  the  reflected  shock  is  not 

5 

exactly  zero.  It  is  several  orders  of  magnitude  less  than  the  ~10 
cm/sec  characteristic  of  Che  fluid  before  it  is  reflected  from  the 
wall.  As  described  above,  these  are  finite  amplitude  sound  waves 
generated  in  the  process  of  reflection  which  decrease  in  amplitude  as 
this  region  expands. 

Once  the  reflected  shock  reaches  the  low  temperature  region 
behind  the  original  contact  discontinuity,  a large  amplitude  weve  will 
propagate  back  to  the  wall  and  modify  tha  uniform  preaaure  and  tempera- 
ture region.  In  the  case  described  above,  we  have  a region  of  low 
temperature  (6o°K) , high  density  materiel  (1,6x10"^  g/cm^)  moving  into 

e region  of  high  temperature  (2300°K) , and  lower  density 
.I4,  5 

(1.0x10  g/cm  ).  For  a long  ahock  tube,  where  the  shocked  region  has 
spread  out  over  a large  distance,  there  can  be  a substantial  time  when 
the  reflected  high  temperature  region  perslate.  Then  except  for  the 
■mall  finlta  amplitude  sound  waves  which  are  present  even  in  experimen- 
tal situations,  tha  medium  can  be  considered  at  rest. 
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At  each  •pedal  and  temporal  location  along  the  length  of  the 
tube  we  can  monitor  the  detailed  production  and  loaa  ratea  for  each 
chemical  apeclea.  For  example,  an  analysis  of  the  chemical 
production  and  loaa  rates  describing  the  55  cm  position  In  Fig.  2 
shows  that  the  chemical  apeclea  are  clearly  out  of  chemical  equili- 
brium. An  extensive  analysis  of  the  rate  scheme  will  be  reserved  for 
a future  paper,  In  which  we  will  simulate  an  actual  laboratory  shock 
tube  experiment. 
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V. 


Conclusion 


In  producing  the  resulte  presented  above,  a premium  has  been  paid 
for  the  complete  time-dependent  prediction  capability.  Both  the 
chemical  klnetica  and  the  shock  hydrodynamics  have  been  treated 
simultaneously  and  accurately  by  explicitly  Integrating  on  the 
appropriate  time  and  space  scales.  Were  we  to  separate  these  calcu- 
lations and  ignore  the  fluid-chemistry  interaction,  we  would  arrive 
at  simplified  representations  which  are  enlightening  but  limited 
quantitatively.  If  we  assumed  ideal  shock  hydrodynamics  and  did 
only  a single,  decoupled  kinetics  calculation  at  one  given  pressure, 
temperature,  and  density,  we  would  not  resolve  the  shock  acceleretion 
due  to  combustion  shown  in  Fig.  2.  Simple  comparisons  such  as  Fig.  4b 
also  illustrate  this  point.  Shock  reflection  off  of  a well  is  another 
situation  where  every  point  in  space  has  a different  chemical  and 
fluid  dynamic  hietory.  The  ideal  conditions  of  sero  fluid  valoclty 
and  constant  temperature  and  prassure  are  altered  by  any  exothermic 
reactions  taking  place. 

We  have  focuaed  our  attention  on  e shock  calculation  because 

it  represents  a severe  test  of  the  resolution  and  accuracy  of  our 

numerical  modal.  This  is  part  of  an  on-going  combustion  program  at 

NRL  to  maasura,  modal,  and  predict  ths  behavior  of  highly  complex, 

interactive  combustion  systems.  We  are  currently  applying  tha 

shock  modal  to  data  from  a shock-initiated  deflagration  which 
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accelarataa  and  lncraasaa  the  shock  prassure  . In  the  future,  tha 
modal  will  be  used  in  coordination  with  shock  tuba  experiments  to 
daclphar  tha  details  of  hydrocarbon  fuel  reaction  rates  end  Induction 
behavior . 


21 


"ti!  ill 


. .:rf  ’W 


i m hi  r ■ 


ACKNOWLEDGMENTS 


This  work  has  been  supported  by  the  Office  of:’  Naval  Research 
through  the  Naval  Research  Laboratory.  The  authors  would 
like  to  acknowledge  the  assistance  of  R.  McCann  and  M. 

Flanigan,  end  the  guidance  and  helpful  suggestions  of  S,  Zalasak, 
E.  Hyman,  D.  Bogan,  R.  Galinas,  T.  Bowman,  D.  Garvin,  and 
G.  Skinner. 


22 


CHEMICAL 


i 


I 

i 


Reacti 

Reactants 

Products 

A 

1 

H 4 H 4 M 

H2  4 M 

1 .60 (-30) 

K 

H + 0 + M 

OH  + M 

2.60(-32) 

B 

H 4 OH 

Hg  4 0 

l.4o(-l4) 

11  + OH  + M 

H?0  4 M 

6.20(-26) 

5 

H 4 H,  0 

fZ 

Hg  4 OH 

1.50C-10) 

6 

H + °2 

0 4 OH 

3.70<-l0) 

i 

H 4 0g4M 

H0?  4 M 

4.10(-33) 

Q 

H + HOg 

H2  +°2 

4.20(-il) 

B 

H 4 H0o 

tl 

0 4 Hg0 

8.3o(-n> 

10 

H + HOg 

OH  4 OH 

4.20<-10) 

11 

H + H;)0£i 

Hg  4 HOg 

2,80(-l2) 

12 

Hg  4 0 

H 4 OH 

3.00(-i4) 

15 

11?  + OH 

H 4 HgO 

3.6o(-ll) 

14 

H2  +°2 

H 4 HOg 

g.io(-ll) 

15 

»2  + °2 

OH  4 OH 

2.80(-ll) 

16 

h2  + ho2 

« + HgOg 

1.20<-12) 

17 

H2  + H02 

OH  4 HgO 

1.21(-12) 

16 

H0  + M 

t— 

H 4 H 4 M 

3.70(-l0) 

19 

0 4 0 4 M 

Og  4 M 

5.00(-35) 

20 

0 + OH 

« + °2 

2.72(-12) 

21 

0 + OH  +M 

HOg  4 M 

2.80(-31) 

23 
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Table  I (continuation) 


WM 

Reaction  Rate  Constants 
k * ATBeC/,T 

Reactants 

Products 

A 

ft 

c 

Reference 

22 

0 + Hp0 

H + H02 

1.75C-12) 

4.50(-0l) 

-2.84(+o4) 

24,21 

83 

0 + H 0 

2 

OH  -4-  OH 

1.10(-10) 

0.00 

-9.24(403) 

19 

24 

0 + HO 

OH  4 02 

fl.OO(-ll) 

0.00 

-5.00(402) 

26 

25 

0 + H2°2 

OH  + HO 

1.40(-12) 

0.00 

-2.13(403) 

21 

26 

OH  + OH 

H 4 H02 

2.00(-ll) 

0.00 

-£,02(+04) 

19 

S7 

OH  + OH 

H2+°2 

1.09M3) 

2. 60 (-01) 

-1.47(404) 

16,24 

28 

OH  + OH 

0 + H?0 

1.00<-11) 

0.00 

-5.50(402) 

19 

29 

OH  + OH  +M 

H202  + M 

2. 50 (-33) 

0.00 

2.55(403) 

19 

50 

OH  + H02 

H£°  + °2 

8.30<-ll) 

0.00 

-5.00(402) 

26 

31 

OH  + H,,0 

\ + H0? 

1.34(-14) 

4.30(-01) 

-3.62(404) 

24,22 

52 

OH  + 0? 

ho2  + 0 

2.20C-U) 

l.ao(-oi) 

-2.82(404) 

24,26 

53 

OH  + HgOg 

HgO  4 H0g 

1.70<-ll) 

0.00 

-9.10(402) 

19 

34 

OH  + M 

H 4 0 4 M 

2.20(-08) 

3.00(-02) 

-5.14(404) 

24,20 

55 

HgO  + 02 

OH  4 HO, 

ci 

2.38C-10) 

1.70  (-01) 

-5.69(+o4) 

24,26 

56 

H?0  + HO,, 

OH  4 H202 

4.70(-H) 

0.00 

-1.65(404) 

19 

37 

H20  + M 

H + OH+M 

5.80C-09) 

0.00 

-5.29(404) 

19 

,58 

o2  + h2o2 

ho2  + ho2 

1.20(-13) 

5.oo(-oi) 

-2.11(404) 

25 

39 

02  + M 

0 + 0 4 M 

3.00(-o6) 

-l.OC 

-5.94(404) 

23 

40 

ho2  + ho2 

°2  + H2°2 

3.oo(-ll) 

0.00 

-5.00(402) 

21 

41 

HO^  + M 

H + 02+  M 

3.50<-09) 

0.00 

-2.50(404) 

19 

42 

H02  + M 

0 + OH  4 M 

1.13(-o4) 

-4.30(-Ol) 

-5.22(404) 

24,25 

43 

H202  +M 

OH  4 OH  + M 

2.00(-07) 

0.00 

-2.29(404) 

19 

* 5 

Rate  Constants  ere  given  in  units  of  era  /see,  T in  °K. 

Three  body  rates  are  in  etn  /jec  vhere  M is  the  third  body. 
Read  1 ,80(-30)  as  1.80  x lo"-'0. 
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Table  IV 
Vlacoaltiaa.  Tfl 


Spaciea 

H 

3.95*10  T* 

«2 

4.71*10”^  T^' 

Me 

7.66x10"^  T* 

0 

l,15xlo“5 

OH 

l.lOxlO'5 

V 

1.60xl0_5 

°2 

1,16x10" 5 

HQg 

1.16x10" 5 

H2°2 

1,1.6x10" 5 

T a temperature  In  °K. 
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Table  V 


Ideal  Shock  Solution 


Ration 

1 

2 

3 

4 

5 

Preciure  (dynae/cn^) 

1. 056x1 07 

1.951*105 

1.951x10s 

1.056x10^ 

Denelty  (g/omS) 

1,076x10"3 

1.557X10"14, 

2.3l4xl0”S 

5.1l8xlo"6 

Velocity  /cm/aec) 

0.0 

1.675*105 

1 ,675x10s 

0.0 

Temperature  (°K) 

300 

g 

6o.3 

1216 

300.0 

Sound  Speed(cm/seo) 

1.0l6xl0S 

1 

4.571xl01*' 

1.009x10s 

5.380X1O1* 

Molecular  Weight 
(amu) 

4 »0 

1 

4.00 

12.0 

12.0 

Gamma 

1.667 

1.667 

i.4o 

l.4o 

Shock  front  velocity  2,05lxl0S 
Contact  dlecontlnulty  velocity  1 . 675*10'’ 
Rarefaction  foot  velocity  1, 218x10-’ 
Rarefaction  front  velocity  -1.016x10'* 
Mach  number  4.00 


cm/eec 


References 


1.  Boris,  J.  P.  and  Book,  D,  L.,  Methods  In  Computational  Physics, 
Vol.  16,  p,  85-129,  Academic  Press,  1976. 

2.  Boris,  J,  P.,  Computer  Physics  Conmuntcations,  12,  67  (I976). 

3«  Boris,  J.  P,,  Flux-Corrected  Transport  Modules  for  Solving 

Generalised  Continuity  Equations,  Naval  Research  Laboratory 
Memorandum  Report  No.  3237,  March  1976. 

4.  Young,  T.  R.  and  Boris,  J.  P.,DNA  I973  Atmospheric  Effects 
Symposium,  Vol.  4,  UNA  Report  No.  31.31  p-9  (VI  75-442),  p.  57I , 
1973. 

5.  Young , T.  R.  and  Boris,  J.  P.,  J.  Phys.  Cham.,  gl,  2424  (1977). 

6.  Getsinger,  R.  W,  and  Schott,  G.  L,,  Shock  Tube  Studies  of  the 
Hydrogen  Oxygen  Reaction  System.  Physical  Chemistry  of  Fast 
Reactions,  Vol.  I Gas  Phase  Reactions  of  Small  Molecules,  B.  P. 
Levitt,  Editor,  Plenum  Press,  1973. 

7.  Bauer,  S.  H.,  Schott,  G,  L.  and  Duff,  R.  E.,  J.  Cham.  Phys.  26, 
1089  (1958). 

6.  Skinner,  Gordon  B.  and  Gorden  H.  Rengrose,  J.  Cham.  Phys.,  43. 

4129,  1965. 

9,  Ripley,  D.  L.,  and  Gardiner,  Jr.,  W.  C.,  J.  Cham.  Phys.  4Jf, 

2285  (1966). 

10.  Gutman,  D.,  Luts,  R.  W.,  Jacobs,  N.  F.,  Hardwidga,  E.  A.,  and 
Schott,  G,  L.,  J.  Cham.  Phys.  4g,  56B9  (I968). 

11.  Galinas,  R.  J.,  Methane  Combust  ton  Kinetics,  Science  Applications, 
Inc,,  Report  SAI/PL-NI-77-02,  I977. 


29 


12.  Stull,  D,  R,  and  Prophet,  H.,  JANAF  Thermochemical  Tablee,  2nd 
edition,  Nat.  Stand.  Ref.  Data  Ser.,  Nat.  Bur.  Stand.,  No.  37, 

June  1971* 

13.  Gordon,  S,  and  McBride,  B.  J,,  Computer  Program  for  Calculation 
of  Complex  Chemical  Equilibrium  Compoeltione.  Rocket  Performance. 
Incident  and  Reflected  Shocke.  and  Chapman- Jouauet  Defonatlone. 
NASA  SP-273,  197 6. 

14.  Hlrachfelder,  J.  0.,  Curtlai,  C.  F.  and  Bird,  R.  B.,  Molecular 
Theory  of  Gaeee  and  Liquids,  John  Wiley  and  Sons,  1954. 

15.  Marrero,  T,  R.  and  Mason,  E.  A.,  Journal  of  Physical  and  Chamlatry 
Reference  Data,  1,,  3 (1972) . 

16.  Williams,  V.  A.,  Combustion  Theory,  p.  419,  Addiaon-Wsaley 
Publishing  Co.,  Inc.,  1965 . 

17.  Bradley,  J.  N.,  Shock  Waves  In  Chemistry  and  Physics,  John 
Wiley  and  Sons,  Inc,,  1962. 

10.  Skinner,  G.  B.,  Mueller,  G.,  Grimm,  V.,  and  Scheller,  K., 
Combustion  and  Flams,  12,  436  (I968). 

19.  Baulch,  D.  L. , Dryadale,  D.  D. , Horne,  D.G.,  and  Lloyd,  A.  C., 
Evaluated  Kinetic  Data  for  High  Temperature  Reactions,  Vol.  1, 
Butterworths,  I972. 

20.  Moretti,  G.,  AIAA  Journal  £,  223  (I965). 

21.  Hampson,  R.  F.,  and  Garvin,  D.,  Chemical  Kinetic  and  Photochemical 
Data  for  Modeling  Atmospheric  Chemistry,  NBS  Technical  Note  066, 
U.S.  Department  of  Commerce/Natlonel  Bureau  of  Standards, 
Washington,  D.C.,  1975* 


30 


22.  Olaon,  D.  B,  and  Gardiner,  W.  C..  An  Evaluation  of  Methane 
Combuatton  Mechanlama.  Paper  preaented  at  the  173rd  National 
Meeting  of  the  American  Chemical  Society,  New  Orleana,  La., 

March,  1977. 

23.  Baulch,  D.  L.,  Dryadale,  D,  D.,  Horne,  D.  G.,  and  Lloyd,  A.  G., 
Evaluated  Kinetic  Data  for  High  Temperature  Reactlona,  Vol,  3, 
Butterworth,  I977. 

24.  Engleman,  V.  S..  Survey  and  Evaluation  of  Kinetic  Data  on 
Reactlona  In  Methane/Air  Combuation.  EPA-600/2- 76-003,  IJ.S.E.P.A. 
Office  of  Reaearch  and  Development,  Waahington,  D.C.,  20460, 
January  1976. 

23,  Bahn,  G.  S.,  Reaction  Rate  Compilation  for  the  H-O-N  Syatem, 
Gordon  and  Breach,  New  York,  I966. 

26.  Lloyd,  A.  C.,  Int.  J.  Chem.  Kinetica  6,  169  <1974). 


31 


TIME  = 0.0  sec 


DIAPHRAGM 

PRESSURE- 

10  atm,  He 

DENSITY  * 

.01  atm,  2H2  ,02  + SMALL  AMOUNTS  OF 

1.076  x lO^g/cm3) 

H,0,0H,H20,H202  ,ho2 

TEMPERATURE' 

5.118  x 10"6(g/cm3) 

300  °K 

300  °K 

J _J 1 

0.0  10.0  20.0  30.0 

DISTANCE  (cm) 

Fig.  1 — Initial  density,  temperature,  and  pressure  for  the 
reactive  shock  calculation 
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TEMPERATURE  <°K) 


DISTANCE  (CM) ► 

Fig.  2a  - Calculated  temperature  ai  a function  of  diftance  after  the  ihock 
haa  traveled  40  cm  down  the  tube.  The  daahed  linei  indicate  reaulte  from 
the  ideal  solution.  The  position  of  the  shock  front  and  contact  discontinu- 
ity are  noted  by  SF  and  CD,  respectively. 
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Calculated  mass  density  as  a function  of  distance  at  the  time 
corresponding  to  Fig.  2a 


TIME  * l.47x  10  SEC 


DISTANCE  (CM) 

Fig.  3 — Profiles  of  number  density  for  H20,  H,  0,  OH,  H20,  and  02  as  a function  of 
position  in  the  shock  tube.  TTiese  correspond  to  the  temperature  profile  in  Figure  2. 


NUMBER  DENSITY  (cm  -3) ► TEMPERATURE  (°K) 


TIME  (sec) ► 

Fit,  4 — (a  I Calculation  of  temperature  change  as  a function  of  time  in  a pre- 
mixed homogeneous  H2-02  mixture  whose  temperature  is  suddenly  raised  to  1100 
(A),  1200  (B).  1500  (C)  CK,  (b)  the  curves  show  water  density  as  a function  of  time 
corresponding  to  the  cuives  in  Fig,  4a,  The  black  dots  are  the  maximum  water 
density  calculated  behind  the  moving  shock  front, 
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TEMPERATURE(°K) 


Fig.  5 - Calculations  of  temperature  as  a function  of  time  at  49  cm.  The 
shock  ic  reflected  at  a rigid  wail  at  50  cm. 
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